Characterization of the channel diffusion in a sodium tetrasilicate glass via 

molecular-dynamics simulations 
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We study the structural and dynamical characteristics of 
the sodium atoms inside and outside the "diffusion channels" 
in glassy Na20-4SiC>2 (NS4) using classical molecular dynam- 
ics. We show that on average neither energetic arguments nor 
local environment considerations can explain the increased 
density of sodium atoms inside the subspace made of the chan- 
nels. Nevertheless we show that at low temperature the mean 
square displacement of the sodium atoms inside this subspace 
is significantly larger than the one of the atoms outside the 
channels. 

PACS numbers: 61.20.Ja, 61.43.Fs, 66.30.Hs 



The diffusion of alkali atoms in silicate glasses is an 
important matter under investigation since several years. 
In particular the properties of sodium atoms inside the 
amorphous tetrahedral network of silica have been the 
topic of both experimental work (!]] and molecular- 
dynamics simulations (MD) since a simple glass like 
Na2 0-4SiC>2 (NS4) can be used as a prototype for more 
complicated glasses. The question of how the alkali 
atoms diffuse inside the glassy network is still a matter 
of debate and a generally accepted theory of ion trans- 
port in glasses is still missing In a previous MD 
study of NS4 we have shown that the ions follow pref- 
erential pathways ("channels") inside the glassy matrix 
H . Nevertheless contrarily to the popular idea proposed 
by Greaves || and developed in further studies [Q] , these 
channels are neither static nor due to a microsegregation 
of the sodium atoms but have to be seen dynamically in 
the sense that the channels are those regions of space in 
which a great number of sodiums have passed during a 
given simulation time. The existence of these channels 
gives rise to a pre-peak in the structure factor at around 
q = 0.95 A -1 seen both in experiments || and classi- 
cal MD simulations H. In this last paper Horbach et 
al. show that the slow dynamics of the sodium atoms 
is closely related to the one of the underlying silica ma- 
trix which is coherent with the fact that there exists a 
strong correlation between the channels and the location 
of the non-bridging oxygen atoms, as shown in previous 
MD studies 

Once the existence of the channels is established, the next 
step is naturally to analyze their characteristics and to 
determine why the sodium atoms take these preferential 



pathways. This is the aim of the present study. In that 
direction we analyze the potential energy and the local 
structure of the sodium atoms whether they are INside 
the channels (Nai„) or OUTside the channels (Naout)- 
It is indeed generally believed that the diffusion of the 
ions occurs via "hopping motions between well-defined 
potential minima" Q which we should be able to detect 
in our simulations. We analyze also the time evolution 
of the sodium densities inside and outside the channels 
which permits us to detect the dynamics of the channels 
as a function of temperature and to quantify the differ- 
ences between the sodium densities, differences that so 
far have only been suggested ||. To elucidate the diffu- 
sion mechanisms of the Na atoms we determine a charac- 
teristic "residence" time inside the channels and we study 
the mean square displacements versus temperature of the 
Nai n and Naout atoms. 

Our system is made of 648 particles, namely Nj$ a = 86 
sodium, 173 silicon and 389 oxygen atoms, confined in a 
cubic simulation box of edge length L = 20.88 A to which 
we apply periodic boundary conditions. The density is 
thus the experimental one, i.e. 2.38 g cm" 3 The 
interaction potential used in our simulations is a mod- 
ified version of the one proposed by Kramer et al |l3|] 
and its complete description can be found in It 
has been shown in previous studies that this potential 
is able to describe reliably many structural and dynami- 
cal properties of different sodium silicate melts ||,[l4| and 
especially NS4 fipT[. We have used as initial structure 
a /3-cristobalite crystal in which we have randomly sub- 
stituted the appropriate number of Si04 tetrahedra by 
Na2 03 "molecules". Subsequently the system is melted 
and equilibrated in the liquid phase at high temperature, 
i.e. 4000 K, for 50000 time steps (= 35 ps). Then it 
is cooled down very rapidly with a linear cooling sched- 
ule at a quench rate of 2.3 x 10 14 K s -1 . During the 
quenching process, the configurations of the system (po- 
sitions and velocities of all the atoms) are saved at dif- 
ferent temperatures (T w 4000, 3100, 2500, 2300, 1900 
and 1700 K). These configurations are subsequently used 
as starting points of production runs of 2 million steps 
(= 2.8 ns) performed in the micro-canonical ensemble 
[(N, V, E) = const]. During these production runs we 
save 2000 configurations equally spaced in time. The 
glass transition temperature of the system can be roughly 
estimated from the location of the bend in a plot of the 
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potential energy vs T to be around T g w 2400 K. At each 
temperature the results have been averaged over three in- 
dependent samples. 

As mentioned above, we have shown in previous stud- 
ies that at low temperature the sodium atoms dif- 
fuse through the (quasi-) frozen silica matrix within a sub- 
set of the total available space. This subset is made 
of pockets (typical size 3-6 A) connected together via 
small pathways (similar in a sense to a neuronal network) 
and this network is what we call channels. Our aim in 
this work is to characterize the channel diffusion of the 
sodium atoms, i.e. to study the possible structural and 
dynamical differences between the Nai n and the Naout 
atoms. In order to determine the channels, we divide 
the simulation box in N to t = 20 3 small distinct cubes of 
volume k 1 A 3 and determine the number density of the 
Na atoms in each of these cubes during the 2.8 ns of the 
simulation. Then we consider only the upper 10 % of 
all the visited cubes. These cubes form "the core" of the 
channels (the space in which the sodium atoms have been 
the more often) or shortly "the channels" . More details 
can be found in Ref. . Once the channels have been 
defined via this time integration, we can go back and de- 
termine for every time t, Nai n and Naout and start our 
investigations. 

Why do the Na atoms visit a specific fraction of the total 
available space more frequently at low temperature? The 
first explanation that comes to mind is that the visited 
sites are more favorable energetically. To check this idea, 
we have calculated the individual potential energy dis- 
tributions of the sodium atoms at 1700 K whether they 
are inside or outside the channels. The two distributions 
have been averaged over the 2000 configurations saved 
at 1700 K and over the 3 independent samples and are 
represented in Fig. 1. It is obvious in Fig. 1 that there 
is no significant difference between the distribution of 
Nai n and Naout^ both are gaussians centered approxi- 
mately at —3.8 eV with a Half Width at Half Maximum 
of w 0.6 eV. It appears therefore that the energy is not 
the driving force explaining why the Na atoms are inside 
the channels. One may argue that the time average can 
smear out the effect but in our previous study we have 
shown that the sodium atoms visit sites previously occu- 
pied by other sodium atoms [|| and therefore if these sites 
were favorable energetically it should appear in Fig. 1. 
We have also calculated the energy distributions of Si and 
O atoms (not shown in Fig. 1 ) and found that they have 
an average energy of —35 eV and —12 eV respectively. 
In fact, the individual potential energy distributions re- 
flect the whole environment of the Na atoms. There- 
fore even if there exists a structural difference at short 
range, it will not be visible in Fig. 1 because of all the 
other (long-range) interactions. It is hence justified to 
look more precisely at the environment of the sodium 
atoms via the radial pair distribution function and the 
integrated number of neighbors. For a given a — f3 pair 
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Note that the number of first nearest (3 neighbors around 
species a, is given by the value of N(r) a -p at the first 
minimum of g a -p(r). In Fig. 2 we have represented g(r) 
and N(r) at T = 1700 K for (a) Na 0u t-Na, Nai„-Na 
and Na— Na pairs, (b) Nao u t— O, Na In — O and Na— O 
pairs and (c) Nao u t — Si, Nai n — Si and Na— Si pairs. The 
major differences can be seen in (a) where it appears 
that a Nai n has on average 1 supplemental Na neighbor 
compared to a Naout atom. This goes together with a 
reduction of the number of nearest Si neighbors (c). On 
the contrary the local oxygen environment (b) seems to 
be the same for the sodiums inside and outside the chan- 
nels. Finally concerning the nearest neighbor distances 
with the other species there is no significant change be- 
tween Nai n and Nao u t except a slight decrease (0.1 A) of 
the Nain— Na distance compared to the distance between 
Naout and the other sodium atoms. It appears therefore 
that the local structure of the sodium atoms inside and 
outside the channels is not significantly different in order 
to explain the existence of the channels. Nevertheless 
it appears from the study of the radial pair distribution 
functions that on average the Nai n atoms have more Na 
neighbors than the Naout atoms. Therefore it is justified 
to calculate at any instant t the number of sodiums inside 
the channels {N^ a (t)) and the number of sodiums out- 
side the channels (N^^ft)). These quantities are shown 
in Fig. 3 at 1700 K (a), 2300 K (b) and 3100 K (c). A 
first information can be obtained from the shape of the 
curves at the different temperatures. At 1700 K (a), the 
fact that the curves become flat after « 0.5 ns indicates 
that the channels are well defined through the frozen sil- 
ica matrix and are made almost of the same small cubes 
during the rest of the simulation. At 2300 K (b), the sit- 
uation is rather different since -/V^j" goes through a maxi- 
mum (and consequently -/V^ "* goes through a minimum) 
after « 1.4 ns. This is a signature of the residual motion 
of the silica matrix which is only guasi-frozen and hence 
a manifestation of the motion of the channels. In fact, 
the channels that we have defined in our calculation at 
that temperature are the common part of the "moving" 
channels wc could have defined over consecutive fractions 
of the whole simulation. It is therefore reasonable to as- 
sume that the real channels at « 1.4 ns (half the simula- 
tion length) are closest to the channels that we have de- 
fined over the whole simulation. At 3100 K (c), the fact 
that the curves are flat indicates that the channels do 
not really exist anymore. What we have defined as chan- 
nels at that temperature are only randomly distributed 
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clusters of density fluctuations. Actually, the matrix is 
melted and the sodium atoms may visit the entire ac- 
cessible space (the system is ergodic) and hence the lo- 
cation of the clusters is completely random and depends 
crucially on the length of the run. Thus, the definition 
of the channels is quite artificial at high temperature. 
Nevertheless when we compare the values of N^ a after 
1.4 ns for the different temperatures we do not see 
an increase of iVj° with increasing temperature contrar- 
ily to what was suggested by Meyer et at || to explain 
an increase of the quasielastic amplitude in the dynamic 
structure factor of sodium disilicate [NS2] melts. Finally 
we see in Fig. 3 that for all the temperatures the major 
part of the sodium atoms is outside the channels. Never- 
theless as shown in || the total volume of the channels is 
small (« 10 % at 2000 K) therefore the sodium densities 
can be very different inside and outside the channels. In 
Fig. 4 we have represented the sodium density pNa in- 
side and outside the channels (which is simply N^ a /V c han 
and A^°" t /(^ 3 " Vchan) where V c han is the volume of the 
channels) as a function of T . It is immediately notice- 
able that the sodium density inside the channels is higher 
than the one outside and in particular at low tempera- 
ture (below T g ). So far we knew that the channels are a 
subspace highly visited by the sodium atoms (they have 
been defined this way) but now we can conclude in addi- 
tion that they are also a subspace of high sodium density. 
We emphasize that a higher sodium concentration does 
not imply necessarily a clustering of alkalis j^J^] . Actu- 
ally we have shown that no clustering can be observed 
on a single snapshot of the simulation box ||. It is also 
noticeable that the two curves converge towards the limit 
of uniform density at high temperature. This is another 
evidence that at high temperature the channels do not 
exist and that what we have defined as channels are ran- 
domly distributed clusters. 

Once we know the position of the sodium atoms with 
respect to the channels it is naturally of interest to cal- 
culate a mean "residence" time of the Na atoms within 
the channels. This can be done by determining the prob- 
ability P(0,t) of a sodium atom inside the channels at 
t = to be inside the channels after a time t. P(0, t) is 
given by: 

, A#k(0) 

P(0,t) = T , - V ndt) (3) 

with: 

J rii = outside the channels 
1 rii = 1 inside the channels 

This probability is represented in Fig. 5 at different tem- 
peratures. Since we save the configurations every 1.4 ps 
(1000 MD steps), the fast decay from 1 at shorter times 
can not be seen in Fig. 5 but since the channels have 



been determined with this "time step" we have decided 
to calculate P(Q, t) in the same conditions. This means 
therefore that our probability is slightly overestimated 
since we do not take into account the vibrational effects. 
In Fig. 5 we see that P(0,t) decreases and converges to- 
wards the long time limit N^/N^ a very rapidly at high 
temperature while at low temperature (T < 2300 K) this 
limit is not reached after 1.5 ns. After subtracting the 
long time limit, these curves can be fitted reasonably 
well by an exponential function with a time constant t 
which is represented versus 1/T in a linear-log plot in 
the inset of Fig. 5. From this inset it is clear that r(T), 
the typical residence time, shows an Arrhenius behav- 
ior with an activation energy around 1.45 eV. This ac- 
tivation energy is larger than the activation energy of 
the Na diffusion process and other more specific charac- 
teristics of the Na displacement (hops between preferen- 
tial sites, forward-backward jumps, decorrelation mecha- 
nism) which is around 1.3 eV jf|. This seems to indicate 
that the mechanism governing the expulsion of the Na 
ions outside the channels is energetically more costly than 
the other mechanisms controlling the mean displacement 
of the ions. 

Finally we have decided to investigate the influence of the 
location of the sodium atoms on their diffusion by calcu- 
lating "conditional" mean squared displacements (MSD) 

(i? 2 (i) = (\r(t) — r(0)| 2 ^}). On the one hand we have cal- 
culated d 2 (t)i D / out which is the "MSD" of the atoms in- 
side/outside the channels at time t = andt ^ (this is a 
reasonable calculation since we know that P(0, t) remains 
non zero over the length of the simulation) . These quan- 
tities are represented in Fig. 6 for (a) 1700 K, (b) 1900 K 
and (c) 2300 K. On the other hand we have calculated 
i? 2 (i)i n /Out strictly speaking: in that case the atoms are 
continuously inside/outside the channels between and 
t. In this latter case the number of steps an Na atom 
is consecutively inside/outside the channels is small es- 
pecially at high temperature therefore this quantity is 
only represented for 1700 K in the inset of Fig. 6 (a). 
From the results represented in Fig. 6 we see that at 
low temperature (T =1700, 1900 K) d 2 n is always larger 
than d 0ut (at long times the two curves collapse since 
the history of the In and Out atoms becomes the same 
on average). This difference is maximum after w 0.1 ns 
when d 2 n is twice as important than c?Qut- This indi- 
cates clearly that the Na atoms are more mobile inside 
the channels. This result is coherent with the findings 
of Horbach et al. |j] who state that the "sodium atoms 
move quickly between preferential sites" but has never 
been quantified directly so far. It is confirmed in the in- 
set of Fig. 6 (a) by the difference between the MSD of 
the Na atoms that are moving only inside the channels 
and the MSD of the sodium atoms moving only outside 
the channels. At T = 2300 K (Fig. 6 (c)), we see that 
the difference between d\ n and d 0ut vanishes even though 
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the channels still exist. Nevertheless as shown in Fig. 3 
at that temperature the channels are not static anymore 
and therefore the distinction between In and Out has no 
clear meaning. At even higher temperatures (not shown) , 
this is of course also true, and the two curves are super- 
imposed. We can hence conclude that when the chan- 
nels are formed and "static" (lowest temperatures), the 
sodium diffusivity is higher when the Na atoms are inside 
the channels. 

In summary, with the use of classical molecular dynam- 
ics simulations on Na20-4SiC>2 systems we have analyzed 
in more detail the properties of the sodium trajectories 
inside the glassy matrix and especially the differences 
between the Na atoms inside the diffusion channels and 
those outside of these channels (the existence of these 
channels defined as the fraction of space mostly visited 
by the Na atoms has been shown in a previous study S). 
Firstly we have shown that the average potential energy 
of the Na atoms is the same whether they are In or Out. 
This is also true concerning their local environment ex- 
cept a slightly higher Na coordination number for the 
sodium atoms inside the channels. This is coherent with 
the fact that, in comparison with the rest of the system, 
the concentration of sodium atoms at low temperature 
(T < 3000 K) is higher inside the channels. At high tem- 
perature the matrix can not be considered as frozen and 
therefore the channels do not exist anymore. We have 
defined a characteristic residence time inside the chan- 
nels t(T) that shows an Arrhenius dependence with an 
activation energy around 1.45 eV. This value is slightly 
higher than the characteristic energy of the other acti- 
vated processes present in this system. Finally, we have 
shown that, at low temperature when the channels can 
be considered as frozen, the mobility of the sodiums in- 
side the channels is larger than the one of the Na atoms 
in the rest of the system: in that sense one can really 
speak about channel diffusion. 

In order to find the origin of this effect and in particular 
the role of the Si02 matrix, one might change the vibra- 
tional characteristics of the matrix or those of the cation 
and see how these changes affect the sodium dynamics. 
These are investigations currently under way. 

We thank W. Kob for interesting discussions. Part 
of the numerical calculations were performed at "Centre 
Informatique National de l'Enseignement Superieur" in 
Montpellier, France. 
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Figure 1 



Figure 3 




FIG. 1. Individual potential energy distributions of the 
sodium atoms inside and outside the channels (1700 K). 
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FIG. 3. Number of sodium atoms inside and outside the 
channels versus time at (a) 1700 K, (b) 2300 K and (c) 3100 K. 
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Figure 5 
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FIG. 5. Probability P(0, t) that a sodium atom is inside 
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of an exponential fit of P(0, t) 



Figure 6 




t[ns] 

FIG. 6. d 2 (i)i n and d 2 (t)out (see text for definition) 
at (a) 1700 K, (b) 1900 K and (c) 2300 K. Inset: Mean square 
displacements R 2 (t)i n and R 2 (t)o ut (see text for definition) at 
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